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Abstract 

Motivated by the ubiquity Dirichlet like energy functionals 
in cellular physics, we study deformations of such systems, 
specifically the constraints governing the interface changes. 
We prove, in general, that the problems of minimisation of 
local mean angular distortion and Dirichlet energy are 
identical (up to a measure). In a more restricted regime of 
hexagonal tessellations, non-existence of minimisers for non- 
linear deformations has been demonstrated. 
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Introduction 



The Dirichlet energy functional, below at (1) is 
common in many areas of physics. 



n 



\Vf\ 2 dx. 



(1) 



Extremal solutions to (1) have been shown to govern 
stable packing's in 2D structures and cylinders 
(Bezdekand Kuperberg, 1990; Cafferelli and Lin, 2007; 
Cybulskiand Holyst, 2005) and implicitly satisfy 
structural stability conditions (Mullins, 1963; Von 
Neumann, 1952). 

Further, in the so-called equal-constant setting in 
liquid crystal theory the Frank-Oseen energy (Oseen, 
1933; de Gennes and Prost, 1993; Virga, 1994) reduces 
to (1). Often, the stability and configuration of such 
structures under deformations/ change of boundary 
conditions are desired. 

This report is concerned with theoretical mathematical 
models to analyse the deformation of 2D cellular 
structures such as one might observe as a cross section 
of tissue, or in various models of foams and elsewhere, 
See (Weaire and Hutzler, 1999; Weaire and Rivier, 
2009). In the latter case, particularly soap foams, the 
theory of conformal mappings has been used to 



various effect (Drenckhan et al., 2004; Mancini and 
Oguey, 2005). In cases where free boundary conditions 
do not always exist, the extension of such approaches 
to fixed boundary conditions, where conformal maps 
do not always exist, is desired. 




FIG 1 EXPERIMENTAL RESULTS FOR BUBBLES TRAPPED 
BETWEEN TWO GLASSPLANES. A. PLATE CONFIGURATION, 
THE BOTTOM PLATE IS HORIZONTAL ANDTHE TOP IS 
SLIGHTLY TILTED. B. COMPUTER GENERATED 
HONEYCOMBSTRUCTURE WITH SHADING PATTERN TO AID 
VISUALISATION OF TRANSFORM EFFECTS. 
C.EXPERIMENTALLY OBSERVED PATTERN FOR THE PLATE 
CONFIGURATION.THE CONFORMAL TRANSFORM FROM 
B -» ClSw = (ia^logiiaz), HERE 27r /a IS THE TRANSLATIONAL 
PERIOD. 

Reproduced from Figure 1 and 3 from (Drenckhan et al., 2004) 

Motivation 

The case of deforming soap films is discussed as it will 
partially motivate the approach we take. A 
deformation of a 2D-cellular structure of soap films is 
achieved in the following manner: consider the 
cellular structure of soap bubbles- finite in extent - 
trapped between two close parallel panes of glass (Fig. 
IB). The steady state with bubbles of uniform size is 
closely approximated by a regular 2D hexagonal 
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tessellation. Vary the angle between the panes. The 
new configuration is expected to be achieved by a 
conformal deformation for two reasons. 





FIG 2 A. SCHEMATIC OF ACCORDION-LIKE HONEYCOMB 
STRUCTURE. B.F-ACTIN LABELLED NEONATAL RAT HEART 
CELLS (CULTURED FOR 1-WEEK). C. MAGNIFIED VIEW OF THE 
CELL CULTURE. THE LABELLING IDENTIFIES STRESS FIBRES; 
THEWHITE ARROWS MARK THE CROSS-STRIATIONS. PD AND 
XD ARE THE PREFERRED AND ORTHOGONAL CROSS- 
PREFERRED MATERIAL DIRECTIONS. THE DIRECTIONS 
ALIGNWITH THE CIRCUMFERENTIAL AND LONGITUDINAL 

FIBRE directions respectively. Scale bars A.200 fxm, B.100 fim,and C. 
10 Reproduced from Figure 2 and 3 of Engelmayr Jr. et al v 2008. 

Firstly, the equations to minimize the interfacial (film) 
energy give, locally, a trivalent equal angle of ^ n /^ 
hexagonal tiling. Secondly, there are no boundary 
constraints. Since the initialconfiguration also exhibits 
/o symmetry, scale invariance andpressure/volume 
considerations suggest that the free boundary "free 
energy minimisation" problem may be solved by an 
angle preserving deformation- that is a conformal 
mapping. Thus it might be expected that the Euler- 
Lagrange equations to be the linear Cauchy-Riemann 
equations. Conformal mappings are also absolute 
minimisers of Dirichlet energy in the free boundary 
case- an elementary consequence of Hadamard's 
inequality. Therefore the appearance of the Dirichlet 
energy functional is expected. Of course the 
connection between conformal mappings and 
harmonic functions is well known. 

Next, in a more constricted regime, one might use 
microfabrication techniques to create accordion-like 
honeycomb microstructure to yield porous, 
elastomeric three-dimensional scaffolds with 
controllable stiffness and anisotropy (Fig. 2) as in 
Engelmayr Jr. et al., 2008. In that paper the authors 
have demonstrated that the neonatal rat cardiac cell 
based tissue, cultured using these honeycomb 
scaffolds, closely matched the mechanical properties of 
rat ventricular myocardium and that heart cell 
contractility was inducible by electric field stimulation 



with directionally dependent electrical excitation 
thresholds showing that honeycombs can overcome 
some of the main structural limitations of previous 
approaches. Here there is a relatively stiff interface - 
scaffolding structure - being deformed at the boundary 
and quite different effects are modelled. 

Objective 

Our aim here is to present a first attempt at a "holistic" 
mathematical formalism which might encompass 
theoretical analysis structures which appear at both 
ends of the spectrum. The film models could be 
extended in various ways. In particular, we should 

• allow the model to include boundary 
conditions on the interface -such as stretching, 
shearing, or other deformations. 

• allow different functionals related to the elastic 
tension of the interface 

• allow for the internal structure of each cell - 
after deformation- to rearrange itself to 
minimise some "free energy" functional(which 
may vary from cell to cell). These functionals 
may range from (1) to fully nonlinear (possibly 
degenerate)variational type energies and could 
reflect possible internal cellular structure. 

Of course, the energy functional characterising tissues 
is rather more complex and we should only expect a 
fairly coarse analogy- however, near steady state most 
symmetries and integral invariants are conserved and 
we will be able to assert existence and regularity of 
minimisers of our models in reasonable physical 
situations. We will also see interesting features 
regarding the singular structure of the solution across 
the interface. 

A discussion on the model's relationship to known 
physics of multicellular structures, the predictability of 
observables such as interface structure, minimal loss of 
symmetry, regularity and so forth under deformation 
will be considered. 

There is one other aspect we will investigate. Many 
naturally occurring functionals have Morrey's 
quasiconvexity property (Morrey Jr., 1966), not 
discussed in detail here, but it is basically an 
assumption that makes the "direct method" in the 
calculus of variations work- minimising sequences 
have a limit which is a minimiser. Roughly it asserts 
the assumption that if a linear mapping is a candidate 
for the minimisation problem (for instance linear 
boundary values), then this linear mapping is the 
minimiser. In 2D, Dirichlet energy obviously has this 
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property as linear mappings are harmonic. 
Quasiconvexityis a natural property of many 
physically motivated functionals, but linear mappings 
are not conformal so some differences are expected to 
be observed. 

The Structure of the Model 

2D-regions in the complex plane tessellated by regular 
hexagons and subjected to a boundary deformation 
are considered. It is assumed that the boundary 
deformation propagates to the interior instantaneously 
to minimise the average angular distortion at 
interfacial vertices. Then the interface relaxes to 
minimise some local energy functional, for instance 
pressure differences across the interface. This then 
moves the interface (and possibly some internal 
cellular organisation) into an extremal configuration. 

Surprisingly we will be able to solve this problem in 
the sense we will identify the initial vertex 
configuration (the average angle minimising 
deformation with given boundary values) and 
discover a curious connection with harmonic 
mappings - though the minimiser will not, in general, 
be harmonic. 

Motivated by the fact that conformal mappings 
preserve angles we consider the case where the 
boundary deformation propagates to the interior to 
minimise the maximum angular distortion. It will 
follow from our analysis below, that this problem has 
a solution (that is a deformation of minimal distortion 
exists) but it may not be unique or particularly regular. 

There are a few things of note in our model. First, we 
should expect the maximum distortion of the 
honeycomb lattice should occur at the boundary. Our 
model will verify this. Second, we should expect that 
our model does not have a completely regular 
minimiser - we should expect some singular structures 
to define the interface. We will also see this in our 
model. 

The Duality Between Local Angular Distortion and 
Dirichlet Energy 

As we indicated in the introduction, conformal 
mappings may arise as energy minimising 
deformations for two reasons. First, they do not distort 
(at an infinitesimal level) angles and second, that they 
are absolute minimisers of Dirichlet energy. 

It is the purpose of this subsection to mathematically 



derive this very interesting relationship. 

In Astala et al, 2010 a connection between distortion of 
relative lengths (linear dilatation) and Dirichlet energy 
was observed. Here we explicitly relate these 
quantities to angular distortion. 

Consider a deformation / of a planar region Q c 
£ .Analytical investigations of deformations require 
some degree of smoothness. Henceforth, it is assumed 
that/ £ M^. 2 (Q,C), that is locally (at each point there 
is a small region where) the first derivatives exist as 
locally square integrable functions. Be aware that such 
a regularity assumption does not even guarantee that 
the deformation is continuous. The reader may 
acquaint themselves fully on the modern aspects of the 
theory of planar deformations in Astala et al, 2010. The 
reader well versed in the analytic theory will see that 
for the most part (and entirely for certain functionals) 
it is possible to work with the much larger, and far less 
regular, class of M^-Sobolev mappings. 

Angular Distortion Functionals 

The local angular distortion is a function of the 
differential map- the best local linear approximation 
to a deformation at a given point. 

To motivate our considerations we start with an 
analysis of the angular deformation by linear maps. A 
linear mapping has the form/(z) = az + bz,a,b £ C. 
The deformation is assumed orientation preserving 
and injective so \a\ > \b\. The injectivity assumption is 
physicallymotivated by the principal of 
"interpenatrability of matter". 

Of course conformal rotation and scaling does not 
change local angular deformation so we may as well 
consider the linear map with \i = a A . 

m2 + fiz (2) 

for our analysis. 

The maximum and minimum angular distortion for 
the linear map are 
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Accordingly, 



max 

0e[O,2n] 

min 

0£[O,2tt] 



= i-IpI- 



(5) 
(6) 



Defining the angular distortion as the ratio of these 
quantities (which would be natural in some sense) has 
a problem for us from the point of view of variational 
calculus. Certainly conformal mappings are absolute 
minimisers as /i = . However [i = is not a 
smoothminimum (in the sense that — ( — ) = 

ds \l—s/ \ s =0 

0, n = 0, and to study variational problems it is almost 
essential that thefunctional has smooth minima so as 
to be able to derive Euler-Lagrange equations and so 
forth. 

For this reason, consider the ratio of 1 + \fi\ 2 and 
1 - I a 1 2 instead. Then, as — f^-V)| = 0, u = will be 

as VI— sV l s= o 

a smooth minimum. 



Since the function t >—> t + 1/t is convex and 
increasing for t > 1 , the maps which pointwise 
minimise the ratio (1 + |(U|)/(1 — |/^|)and those which 
minimise (1 + |/i| 2 )/(l — |/u| 2 ) are identical. 

However, minimising integral averages reveals 
surprising differences. 

Generalised Distortion Junctionals 

The above discussion can be cast in a more general 
framework as follows. 

For a deformation /£ W^ 2 (Q, C) , define the 
Beltramicoefficient 

AO) 



}i/(z) — 



(.'/) 



The theory of quasiconformal mappings and Beltrami 
equations assumes the uniform bounds ||/i(z)|| i «(Q) = 
k < 1 leadingto a "uniformly elliptic" theory. 
However recent advances (Iwaniec and Martin, 2008) 
allow for the consideration of degenerate elliptic 
theory which seems more suitable in situations where 
apriori bounds will not be easy to achieve. 

Note that /i = implies f 2 = (the Cauchy- 
Riemannequations) and the regularity assumptions 
imply / is conformal-angle preserving. On defining 
the distortion measure as the average of the ratios of 
the largest to smallest and smallest to largest angular 
distortion: 
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The calculations applied to the linear differential of 
/show that 

l + l?</(z)| 2 



(9) 



l-|/'/(2)l 2 " 

In fact it turns out that A(z,f) not only controls 
theangular distortion of a mapping, but also its linear 
distortion. 

Example 1. Forz £ Q and < r < dist(z , dCT) ,set 
L(r) = max \f(z) - f(z ) |, t(r) = mm |/(z) - f(z ) \ 

\z\=r \z\=r 

to be the largest and smallest stretching at scale r for 
the mapping/. Define 

L(r) 



K(z Q ,f) =limsup- , . 



(10) 



the linear dilatation of / and, for reasons of 
convexityagain, define a pointwise ratio 



(11) 



It is easy to verify that for a linear mapping z *-* az + 
bzthe number K is simply the ratio of the larger to the 
smaller axesof the ellipse /(S(z ,r)) -the deviation 
from roundnessor eccentricity of the image. 

The following theorem is easy following the 
discussion for linear maps - but nonetheless quite 
remarkable. It shows that the quantities measuring 
angular distortion and deviation from the 
"preservation of roundness" are the same. 

Theorem 1. Let /: £ D — > C be a homeomorphism of 
Sobolev class W^(fl, Q- 

Then 

A(z,/) — K(z,f) almost everywhere in O (12) 
Proof: Using (7), (10) in (11) and simplifying gives 

I/ z | 2 + L£I 2 _I|d/|| 2 



K(z,/) 



A(z,/) (13) 



(8) 



\M 2 -\h\ 2 ' // 

where ||-|| is the Hilbert-Schmidt norm of the 
differentialand Jf is the Jacobian determinant. The 
distortion functionIK is discussed extensively in Astala 
et al, 2009. 

The reader should not be too concerned with the 
"almost everywhere" part of the theorem. It is a 
consequence of working with nondifferentiable 
Sobolev functions - so the quantities involved are only 
defined on a set of full Lebesgue measure. If / is a 
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diffeomorphism, then the equality holds everywhere. 

Eq. (13) provides the first possible connection between 
the Dirichlet energy and angular or spherical 
distortion. This is formalised in the following theorem. 
The proof is a triviality apart from the quite deep 
technical complexity in using the change of variables 
formula for mappings in the stated Sobolev class. This 
change of variables in greater generality than we need 
- for mappings of finite distortion where K(z,f) < 
ocalmost everywhere- has been established by Hencl 
and Koskela, 2006. 

Theorem 2. Let f:f2 — > Hbe a homeomorphismof Sobolev 
class W^' c (n,S2) and supposethat A(z,/) is a locally 
integrable function. Set g=f~ 1 :D. — > D. .Then gE 
W l 1 f{n,n)and 

£a(z,/) dz = JJjDgWW 2 dw (14) 
Proof: Change variables by g in the integral at (14). 



f A(z,/)rfz = / 
JO Jn 

- L 



■ l|D/(z)ll 2 



/(*,/) ^ 
HfJ/t?<™))ll 2 



}(w,g)dw 
J(w,gfdiv 



n J(XMJ) 

f l|D/GfH)|| 2 

Jn J(gW,f)J(™,g) 
fAUf(gW)\\ 2 J&.gfdw 

\\Df{gM)J(w,g)\\ 7 dzv 



- L 

= fj{Dg(u>y)- l J(zi> lS )\\ 2 dw 
Jn 

f \\Dg(zo))\\ 2 dw. 
Jn 



Here, the obvious facts are that f(g(.w)) = w , so 
Df(g)Dg = I andthe Jacobian determinants satisfy 
J (g(w),f)J (w, g) = 1 and the less obvious fact is that 

HD/te(™))/Ks)|| 2 = \\( D 8( W ))~ 1 I( W 'S)\\ 2 



|Wjf(w)| 



(15) 



Note that here it is not the matrices which are the same, 
merely their norm. This fact can be established via an 
eigenvalue argument. 

Observations and More General Models 

Listed below are some observations, theorems, 
assumptions and simplifications that enable the 
development of a general model. 

Existence and Regularity 

Theorem 2 (and in particular (14)) shows that 



minimisers of mean angular distortion have inverses 
which are minimisers of the Dirichlet energy — 
harmonic mappings. In particular we have (following 
§21.1.5 Astala et al, 2010) 

Theorem 3. Let H c M 2 be a convex domain andf 6 

Then the minimisation problem has a unique solution. This 
extremal map is a C K -diffeomorphism whose inverse is 
harmonic inH. 



mm 

ft? 



£a(z,/), f = f <mdCi, (16) 



We sketch a proof of this as follows. First, we may 
assume using the Riemann mapping theorem, that the 
target domain is the disk. This is because if <p: fo(CT) — > 
O (notation explainedmomentarily) is conformal, 
thenA(z,/) = A(z, <p ° /)since <p does not distort angles. 
Therefore we can solve theminimisation problem for 
the boundary values (p ° fo with a smooth 
diffeomorphism / , then the solution we want is 
(p- 1 o / with boundary values f . So we assume/ (Q) = 
D and set g = (fo) -1 . The harmonic extension of 
go\d3 — > 3Q(call it g) is a smooth diffeomorphism 
(given by integration againstthe Poisson kernel), see 
Duren, 2004. Here convexity of the imageD of g is 
needed and the smoothness and non-vanishingof the 
Jacobian determinant is a classical result known as 
Choquet's Theorem (Duren, 2004). Now g minimises 
the Dirichlet energy and so (14) tells us that / = g~ x 
solves our minimisationproblem for mean angular 
distortion. The uniqueness of Dirichlet minimisersfor 
the boundary value problem tells us that / is unique. 

1) Maximum Distortion at the Boundary 

In Martin et al, 2009, it is shown that if /: U — > E 2 is 
a homeomorphism defined in some region and if 
the inverse of /is harmonic, then y. = (i f = f 2 lf z 
satisfies the nonlinear Beltrami equation 

h = i l h ( 17 ) 

Since / will be a diffeomorphism, \fi(z)\ < 1 and 
(17) shows (i to be locally quasiregular (for the 
theory of such mappings, see Astala et al, 2010). In 
particular jx: U — > U 2 is anopen mapping and 
therefore satisfies the maximum principle. That 
is \n(z)\ is largest on the boundary of 11. Since 

the maximum angular distortion also 



A = 

i-IH 2 

occurs on the boundary. 

2) Internal Structures of the Domain 

The above calculations in Section Generalised- 
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distortion-functionals and in particular (14) was 
made for Lesbegue measure. Itcarries over with 
obvious changes in the formula for other 
measures/l (z)dz. That is we have 

jj A(z,/) A(z) dz = jj \\Dg{w)f A(g{w)) dw (18) 

for bounded A. Minimisers of the right-hand side of 
(18)are harmonic with respect to the metric A {z)dz. 
There aremany aspects of this more general theory 
(existence and uniquenessand so forth) that are 
relevant to our investigation in this moregeneral 
setting, see Schoen and Yau, 1997. In particular, 
this mechanism can be applied to relate mappings 
which are harmonic with respect to other metric 
structures and mappings which minimise a 
weighted mean angular distortion. This weight can 
be used to mimic internal structures of the domain. 

3) Other Local Functionals 

In what follows, minimisation of Dirichlet energy 
and harmonic mappings, or elementary 
considerations about pressure/volume and the 
relationship with interfacial curvature will be used 
as the simplest examples of the general models. 
However, there are other natural - and physically 
motivated functionals - which one might consider 
and will lead to a good theory. The prototypical 
such linear functionals are those of divergence 
form: 

£[f] = jj {A(z)Vf,Vf) dz (19) 

with second order Euler-Lagrange equations 
4(z)V/ = 0. Here 4(z) is a measurable positive 
definite symmetric matrix-Dirichlet energy and 
Laplace's equation is obtained for A = I,the identity 
matrix. 

A can be viewed as organising the internal 
structure of a region, as happens in impedance 
tomography and the associated Calderon problem 
(§18 Astala et al, 2010). Here A is viewed as a 
density function and various substructures are 
identified by their relative densities (as a 
measurable function, 4need not vary continuously 
- discontinuous changes in density identify 
interfaces internal to the cellular structure). 

More difficult would be nonlinear examples 
minimising (4(z,/)V/, V/) or related functionals. 
Here, there are significant technical difficulties, but 
some things can be done under stronger 
smoothness assumptions on .4. The linear theory is 
aided by such things as Liouville's Theorem that 



bounded entire, that is defined everywhere on C, 
solutions are constant. And removable singularity 
type theorems, so solutions can be extended as 
solutions over isolated points undermild 
assumptions (boundedness for instance). In the 
linear theory, these things follow from the 
quasiconformal theory by harmonic factorisation 
theorem established in (§16.1.4 Astala et al, 2010) 
which also gives optimal regularity. 

4) Admissible Boundary Deformations 

Our model now describes how the initial 
deformation may propagate into the interior. This 
deformation / with boundary values f gives us a 
"nonlinear" hexagonal tessellation of the interior. 

The procedure will now be used to relax this 
nonlinear tessellation to minimise some free energy 
functional (discussed next). There are two issues to 
consider here. First concerns the local regularity of 
the boundary values. It is a simple consequence of 
a result of Martio and Pavlovic, see (Martio and 
Pavlovic, 2002) and our analysis in section 
Generalised-distortion-functionals relating linear 
and angular distortion, that the angular distortion 
of this extension is bounded if and only if the 
boundary values of /oare bilipschitz. 

That is there must be a constant L so that 

r k - y\ < IA(a-) - fa (.'/) I < L I* - v I for al1 ;/ e an 

(20) 

Roughly this says that the derivative of / is 
bounded and has boundedreciprocal. Strictly one 
needs some smoothness assumptions on dQfor this, 
say smooth and convex. Since there seems little 
physical motivation to consider cases where there 
might be infinite angular distortion we accept this 
condition. Another problem is slightly more subtle. 
Namely a large scale deformation which is locally 
regular may lead to a vertex configuration which 
cannot be relaxed without the interface exiting the 
domain. For instance when/ (Q)is "highly" non- 
convex we expect our relaxation procedure to fail 
as illustrated in Figure 3. 




ABC 



FIG 3 FAILURE OF RELAXATION PROCEDURE FOR A NON- 
CONVEX DEFORMATION. A. DOMAIN Q WITH 3Q COLOURED 
RED. B. A NON-CONVEX DEFORMATION OF Q. C. 
RELAXATION CAUSES EDGES TO SPILL BEYOND THE FIXED 
BOUNDARY. 
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There are basically two ways forward. That is to 
assume the deformation/ has bilipschitz constant 
(the number L in (20)) close to one. Alternatively 
we might assume that the deformation f itself has 
a regular image (say again convex). 

The Interfacial Energy Functional 

The extension of our boundary values to a 
deformation of the domainQgives us a finite mesh 
(collection of vertices and curvilinearedges joining 
them). A connection free energy or interfacial energy can 
be defined as follows. Each vertex v is associatedwith 
a finite number of arcs e h i = 1,2, ■■■ ,n, emanatingfrom 
it (a hexagon model will have n = 3). Let £(e?)denote 
the length of this arc. 

Define 

£l= E EW (2D 

vertices V '=1 

This energy penalises arc length, and it is clear that the 
sum immediately decreases when each arc is replaced 
with the line segment joining its endpoints. 

This is a physically relevant model since the interfacial 
(mean) curvature is a measure of the pressure 
differential across the interface. In an energy 
minimising configuration we expect this pressure 
differential to be so the curvature of the interface 
should be 0, thatis the interface should be locally made 
of line segments (in 2D). 

Here one could explore to useful purpose, other 
increasing functions of the edge length. For instance, 
define O(t) to be increasing on [0, a], O(0) = 0, and 
define the O — connectionfree energy 

by 

E E*(%)) (22) 
vertices <=l 

Again, as O is increasing, arcs can be replaced by line 
segmentsto decrease free energy. 

For a mesh which consists of line segments with 
trivalent vertices, all the internal angles are ^ n /^ at a 
minimum of £ L . 

If not, consider a vertex where this is not the case. 
Choose the smallest of the three angles and move the 
vertex into that region in the direction bisecting that 
angle - keeping all other vertices fixed. A simple first 
order analysis shows that if this angle is not /o, 
then the sum of edge length decreases. A key feature 
here is the linearityof the connection free energy. 



The symmetry of this critical point is not found for 
other functional as at (22), for instance O(t) = t n ,n> 1, 
unlessthere is a very special configuration of edge 
lengths (again easily verified by a first order analysis). 
Further, it is easy to prove that for multivalency, the 
equal angle configuration is a local minimum. 

But for the configuration obtained after a deformation, 
it is not easily seen what the best direction to move a 
vertex is, and it is certainly not always into the region 
with smallest angle (as illustrated below). 




DIFFERENCES FOR MULTIPLE AND 3 VERTEX PROBLEM. 

On the left-hand side, the optimal direction is not into 
the small estangle as the contributions from the many 
other vertices decrease the overall sum. On the right- 
hand side, for three vertices moving into the minimal 
angle lowers energy. 




AN EVAPORATING HEXAGON. 



Another feature of the connection free energy is that it 
allows degeneration. For instance, hexagons can 
evaporate, leaving a vertex with higher multiplicity. 
To see this, consider a regular hexagonal tiling, with 
one of the hexagons, say H, centred at the origin, all 
having unitlength sides. Let 0<A<1 and consider AH. 
The hexagons around H can be extended to six sided 
polygons with all internalangles equal to /o. This 
maintains the local criticality. There are only 12 edges 
affected by this move. Note that the original vertex 
sum is 24 from those edges meeting some vertex of H. 
Replacing H by AH and extending in the way 
described the new vertex sum is 12A, from the edges of 
AH, and a contribution of 12 x (2 — A)from the edges 
that lie in the modified hexagons surrounding AH. The 
vertex sumhas not changed. Though 6 vertices will 
eventually disappear-to be replaced by one. 
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In the same manner irregular hexagons with all 
internal angles ^/gcan shrink, to regular isosceles 
triangles or other polygons all ofwhose internal angles 
are or 71 ® ne can easn y deducea recipe for the 
degeneration (production of defects) as illustrated 
below. 

DEFECTS ARISING. 

One sees that, for instance Pentagons arise from defect 
formation. Of course, if something can disappear 
without energy being gained, then it can also appear 
from certain configurations. 




PENTAGONS APPEAR AND THE INTERFACIAL ENERGY HAS 
NOTINCREASED. 

Of interest is the fact that nonlinear functionals such as 
O(t) = t p ,0 < p and p ^ 1 in (22) do not allow this 
degeneracy and the usual hexagonal packing is rigid 
(no hexagons can be shrunk). 

The sums in question are 12A P and 12 x (2 - X) p which 
are minimised uniquely at X = 1 independent of p 
(aslong as it is not equal to 1). 

There are a few points worth noting. First, if a hexagon 
disappears in this relaxation process, then the 
conformal energy (mean angular distortion) of the 
entire deformation must blow up. Second, if a hexagon 
should appear, then the Dirichlet energy must blow 
up. These two facts are consequences of removable 
singularity theorems for mappings of finite distortion 
with distortion function in some integrability class, 

see Astala et al, 2010 or Marti et al, 2009. For example; 

Theorem 4. Let Q be a planar domain and 11 a compact 
connected subregion which is not a point. There is no non- 
constant mapping of finite distortion f : Q \ 11 — > R 2 which 
shrinks 11 to a point(so f has a continuous extension to ft 
and f (11) =point)with f Q . u A(z,f)dz <<x. 



Of course associated with this observation is the more 
difficult question of precisely what integrability 
classes for the distortion are necessary for the 
conclusion of this theorem, and then the associated 
extremal problems. These are discussed in (Astala et al, 
2010) and the references therein, and so we do not 
pursue the (quite technical) matter here. 

Simulations of Degeneracy 
The data consists of 

- A domain Q tiled by regular hexagons O = 

- A boundary deformation f : dQ — > d£l. 
1) Algorithm 

Step 1. Extend f into the domain Q via thePoisson 
integral of the inverse boundary values so as to 
minimisemean angular distortion. This gives a 
configuration of trivalent verticeslinked by arcs. 

Step 2. Minimise the new configuration by 
evolving the interfacial free energy subject to the 
constraint that the boundary vertices do not move 
(Here, energy minimisation was performed using 
surface evolver (Brakke, 1996). 




A B 




E F 

FIG 4 DEFORMATION SOLUTIONS FOR O(t) = t 11 . A. THE 
INITIALHEXAGONAL MESH, THE BOUNDARY IS COLOURED 
RED. B. ZOOM OF THE CENTREOF D. C. HARMONIC 
DEFORMATIONz >-> z 2 + iz 2 OF A. D. 

2 

RELAXEDCONFIGURATION OF C. E. HARMONIC 
DEFORMATION z >-» z + -z 2 OV A. F. RELAXED 

2 

CONFIGURATION OF E. 
The results of this procedure for some harmonic 
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maps are shown in Fig. 4. At a critical point the 
internal valency is three and the internal angles are 
all near /o and all the interfaces are line 
segments. This gives a "deformed hexagonal" 
packing. 

2) Extending over the Cellular Structure 

At this point we are in the following situation with 
our model. We have 

• A domain D tessellated by regular hexagons 
and a boundary map f . 

• A polygonal mesh in D. of minimum interf acial 
energy, and a correspondence between each 
hexagon H in Q and acell (polygon) P H of 
fi .This correspondencemaps the vertices of 
H £ D to vertices of H 6 U. 

Define a map fn'-dH — > P by the linearextension of 
the vertex values. We now want to fill in the map 
on// by the solution of a variational problem to be 
solved. This problemshould reflect a posited 
internal cellular structure. The filled inmap will 
then be the solution to a second order elliptic 
equationarising as the Euler-Lagrange equations 
for the variational problem. 

The existence of this map can be obtained for the 
sorts of equations described at (19) using the theory 
of c/Z-harmonicfunctions, see Section 6 in Heinonen 
et al, 1993. 

As a first way to advance the model we consider 
the case where there is no internal structure, and 
minimise mean angular distortion (or Dirichlet 
energy of the inverse). In particular we consider the 
case of harmonic mappings. Chouquet's Theorem 
(Duren, 2004) implies that the harmonic extension 
of the piecewise linear boundary values on each 
hexagon is a smooth (real analytic) diffeomorphism 
on each cell. 

The map is not linear in general and the gradient 
does not vanish on the interior, in particular. In fact 
somewhat more is true. First, since each cell is a 
simply connected region we have the 
representation of a harmonic function as 



\H(z)\ = 



<p'(z) 



f{z) 



(26) 



/ H (z) = <p(z) + f{z), (24) 

where <p and ip are holomorphic. The Beltrami 
coefficient,see (7), is 

f'{z) 



}t(z) 



f{z) 



(25) 



but now note that 



Thus the absolute value of the Beltrami coefficient 
is actually the modulus of an analytic function, and 
therefore satisfies the maximum principle. 
Therefore the maximal angular distortion occurs on 
the boundary (it may not in fact be achieved on the 
boundary). This shows that the places where the 
angular distortion occur are confined to the 
interf acial region. 

Theorem 5.The local maxima of the angular distortion of 
the solution f areconfined to the interfacial structure. 

Persistence of the Interface 

The next natural question that arises is whether or not 
the interface is encoded in the global deformation that 
has been found. Basically we ask the following 
question: if we are given a global minimiser of the 
process we have described (filling in the cellular 
structure by harmonic or inverse harmonic mappings), 
is it possible to recover the initial hexagonal packing 
(and therefore its image). We might expect the answer 
is affirmative in general, as for instance the maximal 
distortion is supported on this set. Of course, the 
interface could not be reconstructed if the deformation 
started out with the identity as boundary values. The 
identity would then be the global minimiser and there 
is no structure from which we might recover the 
internal interface. 

However the aim of this section is to show that this is 
about the only case where this occurs. In fact the 
constructed minimiser has the interface as a sort of 
second order Sobolev singularity. 

Theorem 6. Let f be the map deforming a hexagonal 
packing. Then f failsto lie in the Sobolev space W 2,1 across 
each edge of any polygon, unless f is the identity. 

To establish this result we first note that if the interior 
angles are to be preserved, then such a map - if linear - 
will be the identity. Next, the deformation is harmonic 
away from the edges and in the Sobolev space at 
least. Thus Af = away from thevertices by Weyl's 
lemma so in particular / (or it's inverse) will be 
harmonic away from the vertices - so smooth across 
the edges. But the removable singularity theorem for 
harmonic maps (in the easiest form) states that if / is 
harmonic away from an isolated pointand continuous 
over that point, then it is actually harmonic overthe 
point. In particular the singularities cannot be 
confined to vertices. In summary, if the map has 
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regularity across the interface, the mapping is 
harmonic in a neighbourhood of some hexagon. This, 
however, is impossible and will be shown 
momentarily. We want to note here that versions of 
Weyl's theorem is true for a wide class of second order 
elliptic equations (assuming W 1,v - regularityfor some, 
often large, p ) as is the removable singularity 
theorem(which is typically just an integral estimate 
using continuity). 

The following theorem and its proof hold in much 
wider generality, but is used confirm the situation at 
hand. 

Theorem 7. Let T and T be hexagons. Let h:T — > fibe a 
harmonic map with piecewise linear boundary values. Then 
hdoes not extend to a harmonic map of any neighbourhood 
ofPunless h is linear. 

Proof: Suppose there is such an extension. Applying a 
little geometric observation whose proof will leads too 
far astray (but which is clear in the setting) 

Lemma 1. Let Tbe a convex polygon. For v a vertex and 
E v an edge containing v, 

E v = [v + td®, t e [0,s]}, 

with internal angle 8 at v we let Jl v denoteall the rays 
v + te l ® +ke ,k 6 TL.Then there arevertices v x ,v 2 on an 
edge, and rays J2 X £ 3Z Vl and Jl 2 £ R V2 forming a triangle in 
T The situation is illustrated below. 




A TRIANGLE FORMED BY RAYS. 

Returning to the proof, if h has a harmonic extension 
over a vertex, then there is a power series expansion. 
After rotating and scaling (which will certainly not 
affect the statement of the theorem) one may assume 
[0,1] is an edge of both P and ^andft(O) = 0,ft(l) = 1. 
So 

h(z) = az + bz + <p(z), f(z) = £ a k z k + b k z k (27) 

k-2 

A simple first order analysis tells that in fact on the 
two edges with common vertex we must have 
h (z) = az + bz (and consequently b = 1 — a). Thus 
cpmust vanish along these edges. Substitutingy = 
shows 



Q = <p(z)=Y j a k x k + b k x k , 

k-2 

which gives, by uniqueness of power series 
expansions, —a k = S fc and the formula 

<P(z) = I> (z*-Z*). (28) 

k=2 X ' 

Next, for exactly the same reasons as above, <pmust 
vanishon the line {te 1 ®; t £ R}, where O is theinternal 
angle at of J 5 . 

= <p(z) = £ a k t k [e iklt> - e il *J . (29) 

This gives 

a k sin(fc«I>) = 0. (30) 

Observe that the following is immediately implied 

Lemma 2. // the internal angle <f> at any vertex is not a 
rationalmultiple ofn, then h is linear. 

And the proof is simply to observe that in those 
circumstancesa fc = for all k. This was known earlier 
and can be establishedvia repeated application of the 
reflection principle and the observationthat the level 
curves of the harmonic map are real analytic curves. 

Returning to the proof, observe the further implication 
of (28) and (30) that 

h is the same linear map on all of the rays te lk ®, t £ R. 

This of course holds at every vertex. The lemma finds 
an edge and a pair of rays (one might be an edge) 
forming a triangle T on whichft is piecewise linear (but 
linear on each edge). It is then an elementary exercise 
to see that such piecewise linear map always is the 
restriction of a linear map L(z) (not true for general 
polygons, but true for triangles!)- there is a linear map 
from any triple of non-colinear points to any other 
such triple. Then h(z) — L(z) isharmonic on T and 
vanishes on the boundary of the triangle T . A 
harmonic function vanishing on the boundary of a 
simply connecteddomain is identically in that 
domain. Then h{z) — L(z) is real analytic and vanishes 
in the open set T. Thus h{z) = L(z). 

No Global Periodic Minimiser 

In this brief section we consider the large scale 
problem where we have tessellations of the plane 
R 2 and ask about globalminimisers preserving 
hexagonal packings with translational symmetry. In 
the following theorem the groups Tj need not be 
thefull automorphism groups of the packing, just 
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transitive on the cells. 

Theorem 8. Let 7\ and T 2 be periodictessellations of£with 
finite cells and with automorphismgroups T t and T 2 
respectively. So the group sof translations T t permute the 
cells of the tessellation. 

Let C be a fundamental cell in T t . Thereis no solution to the 
minimisation problem for mean angular distortion 

min / A(z,f) dz, (31) 
feJ 7 Jc 

where T consists of all W^ f homeomorphisms /: C — > C 
preserving the tessellations in the sensethat for each 
z h> z + y £ F x there is z >-> z + y £ T 2 such that 

/(z + 7) =fiz)+7> (32) 
unless there is a linear mapping z *-* az + bz + c 
transformingthe tessellations. 

Indeed, there are no critical points at all for the 
variation of (31). 

Proof Since / is a local minimiser, it follows that g is a 
local minimiserfor the Dirichlet functional J||£)(7|| 2 .. 
That is (with theregularity assumptions) g is harmonic. 
Next, the mapping g iseasily seen to have similar 
automorphic properties as / with respectto the groups 
T 1 andr 2 . 

g(z + 7) = g(z) + 7/ 2 z + 7 e Fi, z >-i- z + 7" e T 2 . 

(33) 

Since g is harmonic, g zz = That is g z satisfies the 
Cauchy-Riemann equations and is therefore analytic. 
Notice also that therefore 

gz(z + 7)=gz(z), (34) 

and therefore that \g z \ is periodic and attains its 
maximum value on any cell such asC . Hence g z is a 
boundedentire function which must be constant by the 
classical Liouville'stheorem. The same must also be 
true of the bounded anti-analytic function^. Therefore 
both g z and g z are constanton C and the map g is 
linear. 

Corollary 1. With the hypotheses of Theorem 8. There is no 
solution to the minimisation problem 

jnn JjDf(z)f dz, 

(a harmonic mapping) unless there is a linear mapping 
z h> az + bz + c, a,b,c £ C transforming the 
tessellations. 

Linear Deformations 

Here we examine the case of the linear deformation of 



a regular hexagonalar ray. In models of shearing a 
cellular structure (as indicated below), it is quite 
surprising how much symmetry and regularity is 
retained in the solution. 




FIG 5 DEFORM BY LINEAR MAP, RELAX INTERFACIAL 
ENERGY, EXTEND BY HARMONIC MAP. NOTE THAT WE GAIN 
LOCAL 3-FOLD SYMMETRY 

The procedure Algorithm 1, reduces the sum of the 
internal edge lengths, but the total Dirichlet energy 
has increased (since the linear mapping is harmonic it 
clearly minimises total energy for its boundary values). 

Note. This behaviour is not typical for variational 
minimisation problems. We have already discussed 
how Morrey's notion of quasiconvexity is key and 
would suggest that for shearing the linear boundary 
values should have a linear minimiser. 

One can quite easily see how this symmetry is 
maintained in the shear mapping. As the boundary 
values are linear, the harmonic extension of the 
inverse is also linear. 

As an example, consider the shear {x, y) h> {x + 1, y)on 
the line y = 1 holding the base y = fixed. The 
extension minimisingmean angular distortion (and in 
fact maximal distortion) is the linear map (x,y) ■-» 
{x + y, y) . This moves all the vertical lines on 
whichhexagonal edges and vertices lie to have slope 1 
while preserving the horizontal lines acting as a 
translation on each. 

A local minimum of the functional £/, is found by 
moving vertices upand down lines of slope 1. The 
vertices on the upper boundary remaining fixed. The 
image of an initial cell is then a hexagon all of whose 
interior angles are /o with four opposite edges the 
same lengthand two edges stretched. 

The harmonic extension of the piecewise linear 
boundary values effecting this map between cells 
preserves the slopel effect, by symmetry. That is the 
vertical lines through a pair of vertices are also 
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mapped to lines of slope 1. 

When put together we find the minimiser has the 
property that the all the vertical lines through vertices 
are mapped to lines of slope 1. 

This is a fairly general picture of what happens for 
other linear deformations as well. In order to see 
where the maximal distortion is, and how it is 
distributed, we need to find this harmonic map 
between the cells. This is done in the next section. 

The Mapping on a Hexagon 

In order to better understand the internal cellular 
structure after the deformation, the harmonic maps 
that have been used need to be determined. 
Fortunately in this regular situation this can be done 
explicitly. 

The conformal mapping from the unit disk © = 
{z: \z\ < l}onto a regular hexagon is defined by the 
Schwartz-Christoffel mapping 

1 



JQ 



(1 -Z 6 )V3' 

In this case the image hexagon has side length 



22/3-3 ]-(§) ' 



(35) 



(36) 



and the 6 sixth roots of unity are mapped to the 
vertices of the hexagon. 

The aim is to describe a particular harmonic map from 
ID) to the elongated hexagon T . Suppose the 
boundaryvalues are parameterised by a function 
h :% — > f. 

Then the Poisson integral formula gives 



T(z) 



1 f n 1- 



l-\z\ 



(37) 



From this formula it is evident that if the boundary 
values satisfy the requisite symmetries, then so does 
the conformal mapping. 

Consider a harmonic map from the regular hexagon T 
toT . In order to have the requisite symmetries and 
continuity of the map obtained by extending via theses 
symmetries, the edges cannot be sheared. One can 
expect therefore that the following map is near an 
energy minimiser - it is simply the obvious linear map 
piecewise defined on the edges. 

The hexagon T has all side lengths I and maps to^Pthe 
hexagon with upper and lower side lengths at and 
other side lengths equal to I. The harmonicmap can be 
constructed as follows. 



1) Step 1 

Compute the harmonic mapping defined on the 
disk to T . This is the Poisson integral of 
theboundary values defined by the function h , 



>(e ie ) + | {a -1) 



6 € [-71/3,717/3] 
e' e )-f(n-l) 9 6 [-27r/3,2?r/3] (38) 

c w ) + (it - l)fe($(lf' 9 )) elsewhere 

The map h simply translates the four edges whose 
edge-length is that of the regular hexagon (defined 
by the image of cp), and scales the upper and lower 
edges. 

2) Step 2 

Calculate that (keeping track of orientation) 

<AA%(^ > »)« 

Note that the first term in the last equality is the 
Poisson integral of the boundary values of <pwhich 
is conformal and thereforeharmonic. By the 
uniqueness of harmonic mappings with prescribed 
boundary values, this term must be q> itself 



The next term is in fact the harmonic measure of 
the sectors [— 7r/3,7r/3j U [— n, — 2n/3] U 
[2n/3, n] viewed from the point z. The term is 

rnli 1 - [z[ 2 , „ ( /r 1 - |z| 2 



;r/3 \e" 



T 

JlR 



-dt 



l-|z| 2 



-2 arc Ian 



l-4* + |z| 2 



1~W 



dt 



+2 arc tan 



and sums to 



2 arcfan 



l+ix+\z\ 2 



2y/3(l - \z\ 2 )x 



(40) 



v l-|z[ 2 + |z| 4 -4i 2 , 

Note that one has to be a little careful using these 
formulae in computation because of the branch 
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changes where the denominator vanishes. 
3) Step 3 

Solve the remaining final terms involving 5Re(^>) to 
determine the harmonic map. 

The Harmonic Map 
Starting with the formula 

= f 



dz 



(1_ Z 6)V3- 

Suppose z = e lB with 6 £ [n/3, 2n/3] . Because of 
path invariance of this integral, and if stayed away 
from the sixthroots of unity (where branch changes 
will occur), the following isobtained 



f 

Jo 



dz 



(1-3*) 



6\l/3 



J (1-2^)1/3 + J. J\ 

. r 1 dt r Za 

7 J (l + f6)V3+y. (1 

if (T 



. 2 6)V3 

rfz 



Z 6)V3 



/o (l + f 6 ) 1/3 ' 7; (1 
,. Hi) , f° dz_ 

2V3 r( | } + ( l_ z 6)l/3' 

and so the real part lies in the second term here. 
Assuming 6 > 7r(which can be achieved in general by 
considering the symmetries) and substituting z = e lt , 

ie u dt 



Ji (1 



dz 



(l_ z 6)l/3 



/ 

J 71 



12 (1-^) 



= / 



dt 



dt 



Jk/2 (e 
^ In! 



3it _ g3i'nl/3 

dt 



2 1/3 J nil sm 
cos(30) 

2>1 



1 2 3 



, cos 2 (36) 



2 l/3.3 z 1 ^2' 3' 2' 
where i s ^h e hypergeometric function. This gives 
thefollowing formula. 

Corollary l.For 8 £ [7t/3, 2 7r/3] one /ias 

df 



^ 3 y rc/ 



2 1/3 Jtt/2 sin 
cos(30) 



2i/a . 3 



] ^{3t) 
12 3 



2 3 2 



,cos 2 (39) 



Since at = 5Re(z) = cos(t) and since cos(3t) — 4x 3 — 3x 



<p(z) = - 



x(4^ 2 -3) 
2^3 . 3 



(41) 

for x = x + iy, \z\ = 1, x £ | — i.ijand wherethe sign 
is chosen to be the same as that of y. Actually, fromthis 
and the rotational symmetry of the image hexagon, 



one can deduce the boundary values of the conformal 
map cp at all points. 

The required map is v I / (^~ 1 (z)):r/ — > H. This is 

Yfrt)) =z + g(<p- 1 (z))=z+h(z), (42) 

with h real valued. This shows that the putative 
minimal energymap should preserve the level lines in 
the y-axis. 

Note that 5Re (<p(e i9 )) is the real part of an analytic 
function <p, and is therefore is harmonic. 

Discussion 

Stable configurations of discrete-structures have been 
found under prescribed boundary conditions. 
Assuming that energy stored at every z £ tnt(D), the 
interior to a cell (D) in the discrete structure, is of the 
form 

A(z)Df, and a mappingft £ W 1,2 (d£l, R 2 ), a piecewise 
linear map on the boundary values of the cell and not 
the entire domain, we seek the configuration for which 
the total stored energy (including interfacial terms) is 
minimal. First we find the interfacial configuration, 
then minimise 



/- 

JO 



A(z)\Df\ 2 dz < oo 

over all the class of maps / £ h £ W 1,2 (dQ,]R 2 ) on the 
cell boundary values. For regular v4(z),the existence of 
the minimiser is guaranteed by the principles of 
convex analysis, and the uniqueness depends on the 
convexity of the operator norm. 

It was proved in Theorem 2 that the problems of 
minimisation of local mean angular distortion and 
Dirichlet energy are identical (up to a measure). 
Examples using model energy functionals on 
hexagonal lattices were proposed and their stable 
configurations for linear deformations were 
determined and verified numerically. Note that the 
results, however, are valid for general lattices. 

The results are limited to physical systems that are 
adequately described using mean angular-distortion 
or its weighted-mean-distortion extensions. Note that 
these class of functionals capture a wide range of (non- 
linear) elastic materials. A potential extension would 
be to analyse p — harmonicfunctionals, with Euler- 
equation of the type 

div [\D/r 2 Df) = 0, 

which have been shown to describe nonstandard 
growth conditions (Ruzicka, 2000). 
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There are significant technical difficulties in using the 
direct method of variational calculus to analysing 
these problems for other L p —norms. 

For instance, J Q A p (z,f)dz for p ^ 1 already seems 
completely intractable. This is a subjectfor future 
research. 
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